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Abstract An algorithm for calculating three gauge-invariant helicities (self-, mutual- 
and Berger relative helicity) for a magnetic field specified in a rectangular box is 
described. The algorithm is tested on a well-known force-free model (Low and Lou, 
1990) presented in vector-potential form. 
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1. Introduction 

Gauge-invariant helicities are meaningful measures of nonpotentiality of a magnetic 
field in active regions. These measures in combination with energetic estimates and 
their time dependence may be used for interpretation and prediction of various forms 
of solar activity (see Brown and Priest (1999); Amari et al. (2003a); Amari et al. 
(2003b); Regnier et al. (2005)). 

Self and mutual helicities correspond to the twist and writhe of confined flux 
bundles, and the crossing of field lines in the magnetic configuration respectively 
(Regnier et at, 2005). Time evolution of Berger relative helicity measures the trans- 
port of magnetic helicity through the surface and the effect of boundary transverse 
motions (Berger and Field, 1984). These invariants represent magnetic topology of 
nonpotential magnetic configurations in active regions and its transformation caused 
by eruption. 

This paper considers a problem of exact calculation of the helicity invariants 
within a rectangular box. It is supposed that the full magnetic field vector is pre- 
scribed everywhere in the rectangular box. The helicity calculation problem stated 
in this form as opposed to that stated for a half-space (Longcope and Malanushenko, 
2008) is best suited for fixing physical characteristics of an active region. At present, 
it is possible to calculate force-free spatial magnetic field distribution in a bounded 
volume (see Rudenko and Myshyakov (2009); Rudenko et al. (2010)). Magnetic field 
extrapolation and calculation of helicity invariants and free energy may be used 
for physical analysis of solar activity. In this paper, we examine the algorithm for 
calculating helicity invariants and test it on a force-free model (Low and Lou, 1990) 
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presented in vector-potential form. A key feature of the algorithm is the numerical 
solution to the boundary problem for a potential magnetic field B in the rectangular 
box V in terms of the vector potential satisfying conditions V • A (r) = 0| rG ^ and 
A(r)-S(r) = 0| rGT7 . 



2. Basic formulations and tools description 

2.1. Definitions and general formulation of problems 

Magnetic helicity as a gauge invariant measure may be strictly defined only for 
infinite space or finite volume bounded by magnetic surfaces. The helicity has a 
property to remain constant even when magnetic reconnection dissipates energy. The 
term of magnetic helicity and the condition of constant magnetic topology can not be 
applied to any bounded volume because the vector potential can not be affected by 
an arbitrary additional gradient function (gauge transformation). Choosing potential 
field as a reference field B re f = B pot , 

(Bp t • n)| ay . = (B • n)\ av = g, (1) 



'pot. 



V x A pot , (2) 



one can obtain three gauge-invariant helicity measures (see Berger (1999); Priest 
(1999) and Berger and Field (1984)) for the bounded volume V: 

H seU = / A d ■ B c idv, (3) 
Jv 



H mui = 2 A pot ■ B cl dv, (4) 



v 



A.ff B F = / A cl ■ (B c( + B po4 ) dv- <b C {B c i + B pot ) ■ nds. 



(5) 



Here B c ; = B — B pot ; A c ; = A — A pot ; n is the outer normal to V on S. £ - 
the scalar potential of the gradient field beyond V whose normal component on the 
surface S = dV corresponds to that of A c ; 

(VC-n)ls= (A ci -n)| S! (6) 
where A c ; for B c ; in V is defined as 

Ac (r) = ~ [ x B ci (r') dv'. (7) 

47r Jv |r - r'| J 

Using gauge invariance, one can choose the following conditions on the vector po- 
tential 

V • A pot = 0, (8) 
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(A poi -n)| s = 0. (9) 

Condition (9) is useful for defining reference helicity (5). In this case, the time 
evolution of AHbf describes helicity transport through the surface (see Berger and 
Field (1984)), 

^ ah bf = -2 j> (B A pot ) v pot • nds + 2 j> (v • A pot ) B pot ■ nds. (10) 

The sum of two gauge measures (3) and (4) is the reference helicity defined by 
Finn and Antonsen (1985): 

AH FA {B) = / {A + A pot )- {B- B pot ) dv = H self +H mut . (11) 
Jv 

Reference helicities (5) and (11) coincide only in special cases of V - the half space 
or infinite space above the sphere. 

Given only a magnetic field vector in a bounded volume V, it is necessary to do 
the following steps to obtain helicity invariants (3)-(5): 

a) to solve boundary problems (1), (2) for the vector potential A pot (if A pot is 

known, Eq. (2) gives a solution to B pot for a given g); 

b) to obtain the vector potential of the confined magnetic field A c ; using Eq. (7); 

c) to obtain the scalar potential £ on dV corresponding to (6). 

In this paper, we solve problems a) and c) for the rectangular domain: V = (0, L x ) X 
(0, L y ) x (0, L z ). As mentioned above, problem b) has already its analytical solution 
(7). 

2.2. Problem a) 
Let 

geL 2 (S). (12) 

Consider the problem of finding the vector potential A 6 C°°(V) satisfying equa- 
tions: 

V x V x A = 0, V • A = in V. (13) 
and the boundary condition 

||(VxA°)-n- 5 || i2(g) =0. (14) 

If it is possible to solve the problem (12)-(14), the vector-potential A pot equation 
corresponding to (9) can be written as 

A pot = A — V x A 1 , (15) 

where A 1 is the solution to (13) with boundary conditions 

KVxA^-n- (A°.fi)U| L2(s) =0. (16) 
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The solution to (12)-(14) can be written as a sum: 



A° = ^A 9i +A m . (17) 

i = l 

Here A 9i is the solution to (13) with the boundary function in (14) 

*M=('<v ; r etui tis > 



gds= V x A m -n t ds, (19) 

Si J Si 

where Si - the i-th side of the rectangular box V; A m is the superposition of vector 
potentials given by the system of 5 magnetic monopoles located beyond the box V 

. m \ mtZi x (r — r,-) , . 

A = ? 1 in r 7 — r-^T - 20 

^\T-Ti\[\T-Ti\ + (T-Ti)-Zi] ^ 

This system of monopoles satisfying Equations (19) can be easily constructed by 
fixing their locations and Zj orientations in space and solution of linear algebraic 
system of equations for unknown magnitudes (orts z^ should be chosen so that 
tails of monopoles do not cross the volume V). To each problem (13), (18) on A 9i 
there is a solution giving a unique field V x A 9i , and thus there is a solution to 

(12) -(14) on A , giving a unique B pot . 

Due to the freedom to choose gradient normalization, we can reformulate problem 

(13) , (18) on A 9 ' in terms given by Amari et at (1999): 

AA" =0 in V, (21) 



V • A 9 ' =0 in V, (22) 

Vt ■ Af =0 on S, (23) 

Af = V ± X = V t x x on Si, (24) 

d n A% = on S, (25) 

n, ■ A 9 ' =0 in V, (26) 

-V? X = ^ on Si, (27) 

a„x = on dS 2 . (28) 



Here the subscript t stands for the trace of the operator or the field on the boundary; 
the derivative d n x in Equation (28) means the derivative of the normal to the surface 
dSi in the plane of Sj. 



Paper_hel.tex; 22/09/2010; 0:20; p. 4 



Hclicity invariants of force-free field for a rectangular box 



2.2.1. B VP for the vector-potential function given on one side of the rectangular- 
box 



Further, without loss of generality, is assumed that index i corresponds to the side 
lying in the plane (z = 0). Then A 3 * = from Equation (26). If x is known, tangential 
components A x and A y are BVP solutions to (21)-(26) and satisfy the properties for 
dSi, which follow from Equations (24) and (28): 

A* (x, 0, 0) = d yX (x, 0, 0) = 0, A$ (x, L y ,0) = d yX (x, L y , 0) = 0, 

d x A* (0, y, 0) = d x d yX (0, V, 0) = 0, (29) 
d x A*(L x ,y,0) = d x d yX (L x ,y,0) = 0; 



A* (0, y, 0) = -d xX (0, y, 0) = 0, Af (L x , y, 0) = -d xX (L x , y, 0) = 0, 

d v Af (x, 0, 0) = -d x d yX (x, 0, 0) = 0, (30) 
d y Af(x,L y ,0) = -d x d yX {x,Ly,Q)=Q. 

Properties (29) and (30)) allow us to choose suitable orthonormal bases for A x * and 
A* on S, in L 2 



\J 2 \J T~ cos i 7rmx / L x ) sm("wny/L y ) 



(31) 

0, 1, . . . , oo, n= 1, 2, oo; 



L 77177. 



2-S„ 



— sm(irmx / L x ) cos{nny/L y ), 



m = 1, 2, oo, n = 0, 1, oo. 



(32) 



Here L x , L y are the scales of Si. Basis functions (31) and (32) satisfy (29) and (30), 
therefore expansions A x * and A y * in these bases converge according to the norm of 
W 1 in the neighborhood of dSi (i.e., their expansions converge to them smoothly on 
8Si). 

Using Equations (31) and (32), one can find a harmonical solution in V to problem 
(21)-(28): 



A 9i 



( YZ=o E^i a,™ n h™(x,y) P mn {z) \ 

E"=i E"=i a™h™ n (x,y)p™ n (z) 
V / 



(33) 



Here 



-q m "(2Lz-z) 



sj (Trm/L x ) 2 + (iru/Ly) 2 



(34) 



" d yX {x,y)h™ n (x,y)dxdy (35) 



o Jo 



JO 



d x x(x 7 y)h™(x,y)dxdy (36) 
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Solution (33) satisfies (24) in I? and strictly satisfies (21)-(23), (25) and (26). Choos- 
ing suitable bases (31) and (32) provides smooth solutions in the neighborhood of 
dSi . The magnetic components expressed through partial derivatives A 9i can be 
easily represented by analytical expressions of (33). The same method is applied to 
all sides of V. In numerical implementation of this scheme, expansions (33) may be 
limited by the number of terms corresponding to dimensions of the grid x( x i y) ■ I n 
our implementation, we show integrals of (35) as a sum of analytical integrals in grid 
cells, representing functions d y x and d x \ as a 2D linear interpolation in the cell of 
their grid values. 



2.2.2. BVP-x 



Let us consider the same side (z = 0) as in previous item. The solution to (27), (28) 
will be presented as follows: 

4 

X(r) = x*-X°-J> TJ , (37) 



where 



X* = ^ f g t (r')\n-^—dx'dy' (38) 

is the inhomogeneous solution to (27). The second and third components are poten- 
tial functions satisfying the two-dimensional homogeneous Laplace equation. x° will 
be defined as superposition of the potentials 

X = mix + + m%xy + ni^(x 2 — y 2 ) + mr } (3yx — y )+ 
m 6 (3xy 2 - x 3 ) + m 7 (x 3 y - xy 3 ) + m 8 (6x 2 y 2 - x 4 - y 4 )+ , , 

m 9 (5xy 4 - 10x 3 y 2 + x 5 ) + m 10 {5yx 4 - 10x 2 y 3 + y'')+ { ' 

m u {y 6 - x 6 - 15y 4 x 2 + 15x 4 y 2 ), 

where the vector of coefficients m is selected to satisfy eight conditions of equality 
between tangential derivatives x° an d X* at four vertices of Si and three conditions 
of equality between x° and x* integrals on any three edges of Si. Thus, given that the 
surface integral of §i is 0, the difference of x* —X° have zero tangential derivatives at 
vertcies of Si and zero means on all edges of Si. Let us label Si edges by j and define 
boundary univariate functions Tj equal to values of the normal derivative of x* ~ X° 
on corresponding edges j. x Tj will be defined as potential functions satisfying the 
Neumann problem and the boundary condition: 

\\dnX Tk ~ T k \\ L 2 = on {dSi) k , k = j and . . 

d n = on (dSi) k , k^j. [W> 

Without loss of generality, let j corresponds to an edge on the line y = 0. The 
solution to the problem x Ti m this case can be written in the following form: 



x tj = -e: =1 v*y/£ cobo™* ( °- qm z:^i 



(41) 
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f Lx Pi 
b m = Tj(x,y)t — cos(jtmx/L x )dx. (42) 
Jo V L x 

Solution (41) satisfies (40) in L 2 and strictly satisfies (27), (28) in Si. Solution 
(41) is smooth at the vertices of Si. Components A 9 } and A 9i expressed through 
partial derivatives x Tj can be described by the analytical expressions of (41). The 
same method is applied to every edge of Si. In numerical implementation of this 
scheme, expansions (41) can be limited by the number of terms corresponding to grid 
dimensions Tf c (x,y). In our implementation, we show integrals of (42) as a sum of 
analytical integrals in grid cells, representing functions as a 2D linear interpolation 
in the cell of its grid values. 

2.3. Problem c) 

To obtain potential £ satisfying Equation (6) on the surface S, we use (7). Equation 
(7) yields A c ; on the surface S' limiting the volume V C V 

V = {-Ax,L x + Ax) x {-Ay,L y + Ay) x {-Az, L z + Az), (43) 

where Ax, Ay, Az, are scales of the elementary cell of the discrete grid defined on 
V. To calculate £, we first estimate its values on the grid S' . To do this, we specify 
its zero value at an arbitrary vertex of V Co (r'o) = and then calculate integrals 
on the broken curves I\ jointing vertices r'i on S' and r'o: 

C'(r'i)= f A cl - Ti dl, (44) 
where t is the unit local tangential vector to IV To calculate £ on S, we add integrals 

C(r t ) = C'(r'i) - f ' A c i-ndl. (45) 

J ro 

Integrals (44), (45) are approximated by simple summation. Integrals (44) do not 
depend on Yi curve shape. Therefore in their discrete implementation the choice of 
different configurations of Tj should result in small differences between results of 
integration of (44), (45). 



3. Results of model calculations of helicity invariants 

Let us take the calculation of helicity invariants (3)-(5), (11) as an example for the 
analytical force- free model (Low and Lou, 1990). The expression for the magnetic 
field will be written in terms of vector-potential for n = 1: 

Alz = ~W?+U 1 (46) 
7 r y 1 — A* 

W = -a f -^LdS, (47) 
Jo vl - ^ 
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Table 1. 

H ae ij / H LL H mut /H LL A.H F a/H ll AH F a,ll/Hll AH BFi1 /H ll AH BF2 /H LL 
0.1166 0.6009 0.7175 0.7157 0.7916 0.7810 



where ji = cos 9, 



U=^-, (48) 
r 



P{fJ,) satisfies nonlinear differential equation 

d 2 P 

(l-/i 2 )-7-7 + 2P + 2a 2 P 3 = (49) 

with boundary values 

P(-l) = P(l) = 0. (50) 
Applying operator V x to (46) yields 

Bii=Vx A LL = 1 [ -dgUv - d r U0 + aU 2 $ ) (51) 

that coincides with Equation (3) from Low and Lou (1990). 

For the tests we employed the force-free model calculated in the rectangular box 
V (L x = L y = 1, L z = 0.8; N x = N y = 100, N z = 80) with the following parameters: 
a 2 = a\ l = 0.425; r c = [0.5,0.5,-0.25] are coordinates of the source; the axis of 
dipole lies in the plane (x, z) and has an angle $ = 45° to the 2-axis. 

Next, the notations in (3)-(5) and (11) for helicity invariants are used to denote 
their calculated values. To estimate the accuracy of the algorithm, we have also 
calculated helicity invariant (11) in terms of the known vector potential A^: 



AH FA , LL (B)= f (A LL + A pot ) ■ (B - Bj, ot ) dv. 
Jv 



(52) 



Relative helicity (5) was calculated for two different sets {Ti} (44) and marked with 
indices {1,2}. All helicity invariants were normalized to the true (not invariant!) 
helicity 



H LL = A LL - B LL dv. (53) 
Jv 

Table 1 lists the results. Numbers in Table 1 are helicities of the selected model 
characteristics. The relative errors indicating the level of completion of gauge in- 
variance of Finn & Antonsen reference helicity (3 and 4 columns) and the degree 
of accuracy of Berger reference helicity (5 and 6 columns) are 0.0025 and 0.0134 
respectively. Such sufficiently small errors justify the application of the algorithm 
presented here to calculate helicity invariants of real active regions. 
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